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We study the presence of ferromagnetism in the phase diagram of the two-dimensional honeycomb 
lattice close to half-filling (graphene) as a function of the strength of the Coulomb interaction and 
doping. We show that exchange interactions between Dirac fermions can stabilize a ferromagnetic 
phase at low doping when the coupling is sufficiently large. In clean systems the zero temperature 
phase diagram shows both first order and second order transition lines and two distinct ferromagnetic 
phases: one phase with only one type of carriers (either electrons or holes) and another with two types 
of carriers (electrons and holes) . Using the coherent potential approximation we argue that disorder 
further stabilizes the ferromagnetic phase. This work should estimulate Monte Carlo calculations 
in graphene dealing with the long-range nature of the Coulomb potencial. 

PACS numbers: 81.05.Uw; 71.55.-i; 71.10.-w 



I. INTRODUCTION 

The ferromagnetic instability due to the exchange 
interaction in a three dimensional (3D) electron gas 
attracted attention since the early days of quantum 
mechanics 1 and has been studied in great detail 2 ^. Re- 
cent Monte Carlo calculations 4,5 have confirmed the pres- 
ence of ferromagnetism in the phase diagram of the 3D 
electron gas at low doping. Similar studies have also 
suggested the existence of a ferromagnetic phase in the 
diluted two dimensional (2D) electron gas& with a first 
order transition from a paramagnetic phase to a ferro- 
magnetic phase with full polarization. As the electron 
density is reduced, electron-electron interactions become 
stronger and dynamical screening disappear. At the ex- 
treme limit of zero density the electron gas should crys- 
talize into a Wigner solid where the electrons feel the un- 
screened Coulomb interaction. The elusive ferromagnetic 
phase of the electron gas lurks between the Wigner crys- 
tal and the Fermi liquid state that exists at higher doping 
when electron-electron interactions are fully screened^!. 

In recent years, the experimental search for the ferro- 
magnetic phase of the diluted electron gas has not been 
succesful^SiS. Nevertheless, there has been strong exper- 
imental indications on the existence of ferromagnetism 
in highly disordered graphite samplesiSiii. The origin 
of this phase is still unclear, and a number of different 
mechanisms have been proposecU2ii2*i^*i^. Nevertheless, 
there is no final word on the origin of ferromagnetism 
in graphite. Graphite is a layered material made out of 
graphene layers (a honeycomb lattice with one electron 
per 7T orbital, that is, a half- filled band). The traditional 
view of graphite based on band-structure calculations as- 
sumes coherent hopping between graphene layers, and de- 
scribes graphite as a low density metal with almost com- 
pensated electron and hole pockets, with 1CP 4 to 1CP 5 
electrons per Carbonic. This traditional picture, how- 
ever, completely disregards the strong and unscreened 
interactions between electrons that should exist at low 
densities. In fact, recent experiments in true 2D graphene 
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FIG. 1: Zero temperature phase diagram of a clean graphene 
plane as a function of the coupling constant g, eq. l|T). and 
doping away from half-filling. The dashed line corresponds 
to a first order, and the continuous line a second order 
phase transition between the paramagnetic and ferromagnetic 
phases. The dotted curve corresponds to the value of g with 
eo = 1 and a Dirac-Fermi velocity of hvp — 5.7 eV A, as 
defined by Eq.Q. The points labeled 1-4 in the figure are 
discussed ahead in the text in connection with Fig. |3] 



system^iii&iiiSi show that electron-electron interactions 
and disorder have to be taken into account in order to ob- 
tain a fully consistent picture of grapheneSI. Recent theo- 
retical results^ raise questions on the wisdom of thinking 
of strongly correlated layered system such as graphite, as 
truly 3D. The claim is that the full 2D nature of graphene 
has to be taken into account before graphene planes are 
coupled by weak van der Waals interactions in order to 
form the 3D solid. 

One of the most striking features of the electronic 
structure of perfect graphene planes is the linear rela- 
tionship between the electronic energy, Ek , with the two- 
dimensional momentum, k = (k x ,k y ), that is: e(fe) = 



2 



±fiup|fe|, where vf is the Dirac- Fermi velocity. This sin- 
gular dispersion relation is a direct consequence of the 
honeycomb lattice structure that can be seen as two in- 
terpenetrating triangular sublattices. In ordinary metals 
and semiconductors the electronic energy and momen- 
tum are related quadratically via the so-called effective 
mass, m* , (Ek = h 2 k 2 /(2m*)), that controls much of 
their physical properties. Because of the linear disper- 
sion relation, the effective mass in graphene is zero, lead- 
ing to an unusual electrodynamics. In fact, graphene can 
be described mathematically by the 2D Dirac equation, 
whose elementary excitations are particles and holes (or 
anti-particles), in close analogy with systems in parti- 
cle physics. In a perfect graphene sheet the chemical 
potential crosses the Dirac point and, because of the di- 
mensionality, the electronic density of states vanishes at 
the Fermi energy. The vanishing of the effective mass or 
density of states has profound consequences. It has been 
shown, for instance, that the Coulomb interaction, un- 
like in an ordinary metal, remains unscreened and gives 
rise to an inverse quasi-particle lifetime that increases lin- 
early with energy or temperature^, in contrast with the 
usual metallic Fermi liquid paradigm, where the inverse 
lifetime increases quadratically with energy. 

As mentioned above, its is well known that direct ex- 
change interactions can lead to a ferromagnetic instabil- 
ity in a dilute electron gasi2£. In this work we gener- 
alize the analysis of the exchange instability of the elec- 
tron gas to pure and doped 2D graphene sheets. Al- 
though pure graphene should be a half-filled system, 
we have recently shownSi that extended defects such as 
dislocations, disclinations, edges, and micro-cracks can 
lead to the phenomenon of self-doping where charge is 
transfered to/from defects to the bulk in the presence 
of particle-hole asymmetry. The extended defects are 
unavoidable in graphene because there can be no long- 
range positional Carbon order at finite temperatures in 
2D (the Hohenberg-Mermin- Wagner theorem). Further- 
more, we have also shown that although extended de- 
fects lead to self-doping, they do not change the trans- 
port and electronic properties. Life-time effects are ac- 
tually introduced by localized disorder such as vacancies 
and ad-atoms. Thus, we have also considered the influ- 
ence of disorder in the generation of ferromagnetism. It 
is worth noting that the possibility of other instabilities 
in a graphene plane, related to the Coulomb interaction 
have also been studied in the literature 2425 . The na- 
ture of the exchange instability in a system with many 
bands is also interesting on its own right 2 —, and it has 
not been studied extensively. Furthermore, graphene is 
the basic material for the synthesis of other compounds 
with sp 2 bonding: graphite is obtained by the stacking 
of graphene planes, Carbon nanotubes are synthesized 
by the wrapping of graphene along certain directions, 
and fullerenes " buckyballs" are generated from graphene 
by the creation of topological defects with five and seven 
fold symmetry. Therefore, the understanding of the fer- 
romagnetic instability in graphene can have impact on a 



large class of systems. Finally, we also mention that a 
simple analysis using the standard Stoner criterium for 
ferromagnetism fails in graphene, as the density of states 
of undoped graphene vanishes at the Fermi level2£. 

The electron-electron interaction in graphene can lead 
to other instabilities at low temperatures, in addition to 
the ferromagnetic phase considered here. A local on site 
repulsive term can lead to an antiferromagnetic phase, 
when its value exceeds a critical thresholdSL 2 ^. In the 
following, we will concentrate on the role of the ferromag- 
netic exchange instability, which, as already mentioned, 
is important in electronic systems with a low density of 
carriers, and which has not been considered in the liter- 
ature so far. 

Our main results can be summarized by the zero- 
temperature phase diagram gversus n (where n is the 
doping away from half-filling) shown in Fig. ^ The 
strength of the electron-electron interactions in graphene 
is parameterized by the dimensionless coupling constant, 
g, defined as: 

9 = -r , (1) 

nvF 

where e is the charge of the electron, and €q the dielectric 
constant of the system. Notice that g is exactly the ratio 
between the Coulomb to the kinetic energy of the electron 
system. This coupling constant replaces the well-known 
parameter r s ~ (e 2 /eo)/[h 2 kF /m*} of the non-relativistic 
electron gas (where Uf is the Fermi momentum). In the 
pure compound (n = 0) the paramagnetic-ferromagnetic 
transition is of first order with partial polarization and 
occurs at a critical value of g = g c w 5.3. As the doping 
is increased, the ferromagnetic transition is suppressed 
(a larger value of g c is required) up to around n « 0.2 
where the first order line ends at a tri-critical point a line 
of second order transitions emerges with a fully polarized 
ferromagnetic phase. A unique feature of the ferromag- 
netism in these systems, unlike the ordinary 2D and 3D 
electron gases, is the fact that there are two types of fer- 
romagnetic phases, one that has only one type of carrier 
(either electron or hole) and a second phase with two 
types of carriers (electrons and holes). 

The paper is organized as follows: in the next section 
we present the model for a graphene plane in the contin- 
uum limit taken into account the Dirac fermion spectrum 
and the long-range Coulomb interactions; in Section IlIII 
we discuss the exchange energy for graphene through a 
variational wavefunction calculation in three different sit- 
uations: Dirac fermions without a gap; Dirac fermions 
with a gap; and Dirac fermions with disorder treated 
within the coherent potential approximation (CPA) ap- 
proximation; Section IIVI contains our conclusions. We 
also have included two appendixes with the details of the 
calculations. 
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II. THE MODEL FOR A GRAPHENE LAYER 

The valence and conducting bands in graphene are 
formed by Carbon n orbitals which are arranged in an 
honeycomb lattice (a non-Bravais lattice). The extrema 
of these bands lie at the T point and at the two inequiv- 
alent corners of the hexagonal Brillouin Zone. When the 
filling is close to one electron per Carbon atom, the Fermi 
energy lies close to the corners. Near these points, a stan- 
dard long wavelength expansion gives for the kinetic part 
of the Hamiltonian the expression, 

^w^U-V, kx V kv )> w 

which leads to the dispersion relation, 

e(fe) = ±hv F \k\ . (3) 

In a tight-binding description of the graphene plane with 
nearest neighbor hopping energy t the Dirac-Fermi ve- 
locity is given by: 

3 

hv-p = -ta (4) 

where a is the Carbon-Carbon distance (t ss 2.5 eV and 
a = 1.42 A)i£. The eigenstates of (01 can be written as: 

*w(r) = ( v ;; ( ( r r j)x., 

_ ^ / e # fe /2 X 

I 



where a and b label the two sublattices of the honeycomb 
lattice, (f>k — arctan(fc a /fc x ) is a phase factor, a = ±1 
labels the electron and hole-like bands, and x<r is the 
spin part of the wavefunction. The dispersion and the 
wavefunctions are the solutions of the 2D Dirac equation. 
This approach in the continuum requires the introduction 
of a cut-off in momentum space, k c , in such a way that 
all momenta, k, are defined such that: < \k\ < k c , 
where k c is chosen so as to keep the number of states in 
the Brillouin zone is fixed, that is, nk 2 = (27r) 2 /A , and 
Ao is the area of the unit cell in the honeycomb lattice. 

It is easy to show that with the dispersion given in 
© the single particle density of states, p(E), vanishes 
linearly with energy at the Dirac point, p(E) oc \E\. In 
this case, there is no electronic screening^ and the elec- 
trons interact through long-range Coulomb forces. The 
electron-electron interactions can be written in terms of 
the field operators, ^(r), as: 

Ul = \\ dr i dr 2* t ( r i)* t (^2)F(n-r 2 )*(r 2 ) 1 i'(r 1 ), 

(6) 

where V(r) = e 2 /(eor) is the bare Coulomb interaction. 
One can now expand the field operators in the basis of 
states given in (|SJ), that is, 

*(»0 = -4= E ^k, a A r ) a k, a ,<j (7) 

V A . 

k,a,{T 

where ak, a ,a (a^ „ a ) is the annihilation (creation) oper- 
ator for an electron with momentum k, band a, and spin 
fj (a =|, I and A is the area of the system). In this case, 
the Coulomb interaction reads: 



k,p,q cti ,. . . ,c*4 (T,cr r 



(8) 



It is easy to see that the Coulomb interaction induces 
scattering between bands (inter-band) and also within 
each band (intra-band) . Furthermore, the 1/q depen- 
dence of the interaction (that comes from the Fourier 
transform of the 1/r potential in 2D) provides an 
electron-electron scattering that is stronger than in 3D, 
allowing for the possibility of a ferromagnetic transition 
at weaker coupling. As in the case of the Hund's coupling 
in atomic systems, the spin polarized state is always pre- 
ferred when long-range interactions are present since, by 
the Pauli's exclusion principle, both kinetic and Coulomb 
energies are minimized simultaneously. This should be 
contrast with ultra-short range interactions of the Hub- 
bard type that almost always benefit anti-ferromagnetic 
coupling via a kinetic exchange mechanism. 



III. EXCHANGE ENERGY OF A GRAPHENE 
PLANE. 



In what follows we examine the required conditions 
for a ferromagnetic ground state in graphene. Our pur- 
pose in this work is not to obtain exact values for the 
critical couplings, that may required more sophisticated 
approaches, but instead our aim is to show that a ferro- 
magnetic ground state in graphene is possible in princi- 
ple. In order to study the ferromagnetic instability we 
use a variational procedure that respects all the sym- 
metries of the problem. We assume that: (i) the ferro- 
magnetic instability only affects states close to the Dirac 
points in the region at the edge of the Brillouin zone (that 
is, long wavelength approximation is still valid); (ii) in 
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the ferromagnetic state the electronic bands are shifted 
rigidly (hence, self-energy effects such as Dirac-Fermi ve- 
locity renormalizations are neglected); (Hi) even when 
the bands are shifted, and a finite density of states is 
produced at the Fermi energy, the Coulomb interaction 
remains unscreened (this assumption is equivalent to as- 
sume that the chemical potential shift is always small 
and that the screening length is larger than the inter- 
particle distance); (iv) the ferromagnetic state is uniform 
and translational invariant. Besides considering the case 
of a gapless system, we have also studied the case where 
a gap A opens in the Dirac spectrum (that is, when the 
dispersion relation becomes = ±yA + h 2 v F ). The 
gapped case is interesting because it allows the study of 
the crossover between the Dirac case when A = to the 
standard 2D case with a finite effective mass m* oc A 
(see details ahead). We also briefly the discuss the ef- 
fects of disorder on the stabilization of the ferromagnetic 
state via a CPA approximation in order to point out that 
disorder may be fundamental for the realization of a fer- 
romagnetic phase in graphite. 



A. Gapless system 

1. Exchange energy. Inter- and intraband contributions. 

The possible ferromagnetic instability arises from the 
gain in exchange energy when the system is polarized. A 
finite spin polarization, on the other hand, leads to an 
increase in kinetic energy. Thus, there are two compet- 
ing energies in the problem: the exchange energy that is 
minimized by polarization and the kinetic energy that is 
increased by it. The variational states that we consider 
in our approach are Slater determinants of the wave- 
functions given by (J3J) in the configurations shown in 
FigH 

As function of the Fermi wave vector, k-p, the kinetic 
energy of the unpolarized state is: 
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and the exchange energy, for any doping, as determined 
from Eq.(jHJl can be written as 
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where rip Q ° (k) is the Fermi occupation function, a(b) 
is the band indice, and a a , otb = ±1. 

In the ferromagnetic state the degeneracy of the spin 
states is lifted and the Fermi momentum of the up and 
down spin states becomes fcf and k\, respectively. De- 
pending on the values of k-p^k^ and fci, we can define 
the three cases shown in Fig. For a doping, S per unit 
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FIG. 2: Occupied and empty states in the paramagnetic and 
ferromagnetic ground states of Dirac fermions (a) half-filling 
case; (b) finite doping and one type of carrier in the ferro- 
magnetic phase; (c) finite doping and two types of carriers in 
the ferromagnetic phase. 



area, the number of electrons per Carbon away from half- 
filling, n, can be written as: 



SA 



11) 



Because of the different values of fcj and k^ the system 
acquires a spin magnetization, /i = g s ^Bm, where g s ss 2 
is the electron gyromagnetic factor, (Ib is the Bohr mag- 
neton, and to = s Ao with s = n-f — nj,, is the spin polar- 
ization. Notice that the maximum polarization allowed is 
777 = 2 — 2n since each added (subtracted) electron leads 
to a doubly (empty) Carbon 7r orbital. 

The total exchange energy, ea. Hl(J|) . can be split into 
intra- and inter-band contributions. In many band sys- 
tems where the different bands arise from different atomic 
orbitals, the overlap integral between Bloch states corre- 
sponding to different bands can be neglected, and, con- 
sequently, there are no inter-band contributions to the 
exchange energy. An analogous effect arises when the 
different bands are localized at different sites of the lat- 
tice, as in the gapful case to be considered below. There 
are also situations where the different bands arise from 
the same orbitals at the same sites, but their phases in 
a region much larger than the unit cell are such that the 
overlap integral vanishes. This is the case for the two 
different Dirac cones which can be defined in the honey- 
comb lattice. We do not need to include in eq. I|1(J|I terms 



5 



due to interactions between electrons near different Dirac 
points of the Brillouin Zone. 

The case studied here, where the overlap between 
Bloch states in different bands cannot be neglected, and 
a corresponding term in the exchange energy has to the 
included is generic to narrow gap semiconductors, and 
this term may be important in lightly doped materials. 

It is worth noting that these inter-band exchange effect 
arise from the non local nature of the exchange interac- 
tion. They cannot be studied when the exchange energy 
is approximated by a local term which only depends on 
the total charge density. 



2. Undoped case: n = 



phase, and fcj = ^2n(s + S) and k± = ^2n(s — 5) in 
the ferromagnetic phase. In this phase there is only one 
type of carriers, either electrons or holes. The change 
in energy between the paramagnetic and ferromagnetic 
phase is: 

AE = AK + AE CX = —v F h(k? + fcf - 2fc|) 

(Z,TT ) €q 



2k c klR 2 



4 k c k F R 2 
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kc 
k-p 



2k c k 2 R 2 ' Kl 



(14) 



The Fermi level in the paramagnetic case is at e F = 0, 
and the bands are half-filled. Then, in the paramagnetic 
state one has fcj = fci . When the system polarizes the 



magnetization is such that fcf = \Z2tts and the change in 
energy relative to the paramagnetic state is given by: 



AE 
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(2tt) 2 4e 
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4fc c A;?i?o 



,(12) 



where the functions R n (x) are defined in the Appendix 
lAl Unfortunately it is not possible to find an analyti- 
cal expression (using elementary functions) for the en- 
ergy change as a function of the electron polarization 
s = fc|/(27r). For fcf <C k c , the leading contribution 
comes from the expansion of function Rq(x) ~ — xln(x) 
for x <C 1 (see Appendix 0. Hence, the exchange en- 
ergy increases as the polarization increases, and a ferro- 
magnetic state with small magnetization is not favored. 
This effect can be cast as a logarithmic renormalization 
of the Fermi energy, which reduces the density of states 
near the Fermi level, and suppresses the tendency toward 
ferromagnetisms*. 

At large magnetizations, k 2 /s ~ 1, the kinetic energy 
contribution tends to a term proportional to v F k^ and 
the exchange contribution becomes negative and propor- 
tional to — (e 2 kl)/eo. The exchange term dominates, 
and the system undergoes a discontinuous transition to 
a state with polarization of order unity when: 
r ;2 



9c 



> 



16tt 



hv F e ~ 6i?i(l) - 12i? (l) 



5.3, (13) 



5.3 for 



which gives the critical coupling g c {n = 0) 
the appearance of ferromagnetism in the clean system, 
as shown in FigQ] 



3. Doped case, one type of carrier in the 

ferromagnetic phase 



The behavior of the energy change as a function of the 
spin polarization is shown in left hand pannel in Fig[3| 
for points 1 and 2 of the phase diagram in Fig^ No- 
tice that the transition between the paramagnetic phase 
(point 2) to the ferromagnetic phase (point 1) is discon- 
tinuous with full polarization, m = 2 — 2n. In this case 
analytical expansion when s <C 5 is now possible. For 
&f > fc-f , k[ -C fc c the value of the exchange contribution 
is dominated by the expansion of R 2 (x) (see Appendix 
|A"}| . The contribution of the exchange interaction to the 
term proportional s 2 is positive at low doping, and a 
continuous ferromagnetic transition is not possible. This 
contribution becomes negative only for n — SAq > 0.059. 
As in the previous case, we can also analyze the system 
energy for large values of the magnetization. We obtain 
an instability to a ferromagnetic state with full polar- 
ization (rn = 2 — 2n), which for « / leads to a state 
with both electron and hole carriers with different Fermi 
surface areas. The dependence of the coupling constant 
g c on n is given in Fig^ by the dashed line. (See more 
on the conclusions about a speculative scenario for the 
origin of electrons and hole pockets in graphite.) 



Jf.. Doped case, n^O, two types of carriers in the 
ferromagnetic phase 

In this case the calculation is analogous to the previous 
one. The change in energy in this case is given by: 

AE = AK + AE CX = — v F h(k? - - 2kl) 

D7T 1 i 



A e 2 
(2^) 2 eo 

2fc c fc 2 i?i 



***** I T- 



[fcf.Ri(l) + fcji?i(l) - 2fcf-Ri(l) 



2k c k 2 R 2 1 Kl 



(15) 



In this case the doping, 6, and magnetization, s, are 
such that &f = V2nS in the paramagnetic paramagnetic 



As in the two previous cases, the leading term when 
fcp,fc-|-, fej <C k c is due to the expansion of the function 
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i?2(x), which leads to an increase in the exchange en- 
ergy, which is detrimental for ferromagnetism. The en- 
ergy change as a function of m is shown in the right hand 
panel of Fig|3J We show the energy at points 3 (param- 
agnetic) and 4 (ferromagnetic) of Fig. ^ The transition 
in this case is second order with only partial polarization, 
m < n. As a consequence only one type of carries exist. 
The dependence of the coupling constant g c on n is given 
in Fig^Jby the solid line. 




0.5 1 1.5 

spin polarization m 



0.5 1 1.5 
spin polarization m 



FIG. 3: Behavior of the energy curves as function of the 
magnetization for the points marked in the phase diagram of 
Figd 



B. Gapful system 

A gap can open in the Dirac spectrum when the two 
sites in the unit cell of the honeycomb lattice model be- 
come inequivalent equivalent. In this case, the kinetic 
energy Hamiltonian, Eq.(J2J changes to: 



Wkin(fc) 



A 



v F h(k x + iky) 
-A 



which leads to the modified dispersion relation, 



e fc = ±v/A 2 + (^ F |fc| 



(16) 



(17) 



For wavevectors such that /wF|fc| 3> A the energies and 
wavefunctions are essentially the ones found in the ab- 
sence of the gap, as discussed previously. If the filling is 
such that the Fermi wavevector satisfies this conditions, 
but kp <C k c the analysis presented earlier remains valid. 
At sufficiently low fillings, Ti^fI^f! ^ A, the dispersion 
relation, Ea. (|16|l can be approximated by: 



±A± 



(^ F |fc|) s 
2A 



(18) 



and the bands depend quadratically on the wave vector 
and we can define an effective mass m* = A/v F . Hence, 



the contribution of the kinetic energy to the polarization 
energy is formally similar to that obtained for an 2D elec- 
tron gas with parabolic dispersion discussed extensively 
in the literature. In this case, the spinor wave function 
becomes: 



*fe,<r( r ) 







Xo 



(19) 



for the upper sub-band, while the weight of the spinor 
is concentrated on ij) , Eq.lp]!. for the lower sub-band. 
This change modifies significantly the spinor overlap fac- 
tor in the calculation of the exchange integral, Ea. (|10fl . 
The overlap between Bloch states in different bands for 
momenta near the Fermi points vanishes (see the discus- 
sion at the end of Section III.A.l). These states do not 
give rise to inter-band contributions. The only inter-band 
contributions which need to be included arc due to inter- 
actions between states far from the chemical potential 
among themselves, and between these states at the bot- 
tom of the lower band and those at the Fermi level. These 
terms are not modified when the system is polarized, and 
they do not contribute to the exchange instability. The 
remaining intraband term is equivalent to that derived 
for the electron gas with parabolic dispersion relation. 
The change in energy when the polarized state is formed 
can be written as 



AE 



AK + AE r 



Stt 2A [k 



kf - 2fc F ) 



A 



(2tt) 2 e 3 



=(*? + *; 



2*|). 



(20) 



As in the usual case of the 2D electron gas, the system 
shows an instability toward a ferromagnetic state when 
fcp < (16Ae 2 )/(7TUp£o)- In agreement with the previous 
discussion, this instability vanishes when A — * 0. 



C. The effect of disorder 

We approximate the effects of disorder on the average 
electronic structure by means of the CPA 31 . This approx- 
imation describes the effects of disorder on the electronic 
structure by means of a local self energy, which is 

calculated self consistently. While CPA cannot describe 
localization effects, it still gives very good results for the 
physical properties of graphene^i. 

The total energy, including the exchange contribution, 
can be expressed in terms of single particle Green's func- 
tions, which are calculated within the CPA. The main 
steps of the calculation are sketched in Appendix B. Wc 
assume that the disorder is induced by vacancies, as likely 
to occur in samples treated by proton bombardment. The 
amount of disorder is parametrized by the concentration 

vacancies, n vac . The CPA leads to a density of states 
which is finite at u> = 0, and decays for u> 3> i>F?iva<~. 

Assuming that lim^^o ImS(w) = So ~ (hv-p)/l, where 

1 is the average distance between vacancies^ the calcu- 
lations in Appendix B admit some simplifications. If 
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the concentration of vacancies is small, Eo <C hv-p\k c \. 
At large energies the CPA result vanishes quite fast as 
a function of energy, lima,— >±Su F |fc c | = 0. Disorder 

only changes significantly the results obtained for a clean 
plane if £f *C So- This regime corresponds to electronic 
densities such that \n\ <C n = {H /hvp) 2 /2ir. 
In this limit, we can approximately write 



fi«F|fc|>So, 

1/2+^ ^ F |fc|«E , 



(21) 



where the ± index refers to the two subbands of the non- 
interacting system (see Appendix B). 

The total density of carriers is obtained by integrating 
this expression over k (see Appendix B). Finally, we can 
also calculate the density of states per unit area and unit 
energy, which, for < Eq, becomes a constant: 



D{u) = D a 



1 So, / 



2tt vl 



|cf| < S , 
(22) 



A constant density of states implies that the total number 
of carriers scales as n ss D ep, instead of the relation 
n oc e| obtained for the clean system. 

From equations (|21|) and l|22|) we can infer that both 
the kinetic energy and the exchange energy depend 
quadratically on the density of carriers, since K(n)—K(0) 



and E cxch (n) - E cxch (0) scale as e|(n) 



In addition, 



we know that for n ss no the values of K(n) and -E xch(?i) 
should be comparable to those obtained in the absence 
of disorder. Then, we can write: 



K{n) 
E exc h(n) 



Ckii 



2A ^ 
3irh 2 Vp \no 



Ccxch 



A e 



2y3 
^0 



"0 



(23) 



where Ck; n and c cxc h are numerical constants of order 
unity. In a spin polarized system, we have: 



E tot (n, m) = - [K(n + m) + K(n - m) + E cxch (n + m) + E cxch (n - m)] 



so that: 



AE = AK + AE cxch = c kin 37r ^2 I — ) ''-nru 



2A T.l ( m 
I — 



(24) 



(25) 



The ferromagnetic phase is stable provided that: 

27TC kin 



fjc, disorder 



> 



Ccxch 



(26) 



This result implies that, if n <C riQ the critical coupling 
is independent of the amount of disorder. 

We have estimated the ratio c cxc h/cki n performing nu- 
merically the calculation described in Appendix B for 
suficiently low carrier concentration and density of va- 
cancies. We find: 



^?c,disordc 



3.8. 



(27) 



indicating that in the case of disorder ferromagnetism 
is stabilized at a smaller value of the Coulomb interac- 
tion. Thus, we can conclude that, at least in CPA, fer- 
romagnetism will be enhanced when disorder is present, 
in agreement with the experimental data— 

The enhancement of the tendency towards ferromag- 
netism in the presence of disorder is due to the increase 
in the density of states at low energies. The existence 
of these states implies that a finite polarization can be 
achieved with a smaller cost in kinetic energy, in a qual- 



itatively similar way to the Stoner criterium which ex- 
plains itinerant ferromagnetism in the presence of short 
range interactions. 



IV. DISCUSSION AND CONCLUSIONS 

We have analyzed the ferromagnetic instabilities in- 
duced by the exchange interaction in a system where the 
electronic structure can be approximated by the 2D Dirac 
equation, as it is the case for isolated graphene planes. 

In pure graphene we have found that, as a function of 
doping, a ferromagnetic transition is possible when the 
coupling constant is sufficiently large. Our findings are 
summarized in the zero temperature phase diagram pre- 
sented in Fig. n I n this figure we represent the critical 
coupling g c as function of the doping n. There are two 
different regions in the phase diagram. For small doping, 
n < 0.2 the transition is first order, leading to a ferromag- 
netic phase with spin polarization m = 2 — 2n and two 
types of carriers (electrons and holes). For doping larger 
than n > 0.2 the transition becomes of second order with 
a magnetization smaller than the doping n and one type 
of carrier (electrons or holes). The connection between 
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the magnetization and the carrier type is unique to the 
Dirac fermion problem. We should emphasize that our 
calculation for the Dirac fermion problem is at the same 
level of the one performed by Bloch, and therefore it is to 
be expected that an exact solution of this problem will 
modify quantitatively the phase diagram analyzed here. 
It is also worth remarking that the electronic structure 
shown in panel of (c) Fig. [SJshows that, in the ferromag- 
netic phase, a nominally half filled system has electron 
and hole pockets. The existence of these pockets does not 
depend on the presence of intarlayer coherence, however 

We have also analyzed the effect of the exchange inter- 
action in disordered systems using the CPA. A continu- 
ous transition into a ferromagnetic phase is possible, and 
the coupling required for its existence is reduced with 
respect to the clean case. This tendency can be qualita- 
tively explained by noting that the disorder leads to an 
increase of the density of states at low energies, making 
the system more polarizable. This explanation is rather 
general, and it should not depend on the way the effects 
of disorder are approximated. 

Finally, one would ask how our results can be trans- 
lated for the experiments in disordered graphit e) 10 ' 11 . If 
we naively think of graphite as a stacking of isolated 
graphene planes we can estimate the value of the coupling 
constant for graphite to be g ~ 2.8 (for eo rs l) 1 ^, and 
therefore far away from the ferromagnetic region (corre- 
sponding to the dotted line in FigCJ. The presence of 
disorder will definitely bring the value of the critical cou- 
pling to lower values and according to our calculations 
5c, disorder ~ 3.8 would put dirty graphite at the border- 
line of a ferromagnetic instability. 

Nevertheless, the picture of graphene as a non- 
interacting stacking of graphene planes is certainly in- 
correct. Because of the absence of screening, long-range 
forces will play a major role, and the graphene planes will 
interact via van der Waals interactions. The problem of 
ferromagnetism in graphite still depends on the better 
understanding of the coupling between graphene planes. 
More work has to be developed in order to understand 
the problem of ferromagnetism in graphite. In any case, 



our results here are valid for single graphene planes and it 
would be very interesting to investigate whether graphitic 
devicesiii£ii2i2& studied recently can sustain any form of 
ferromagnetism. 
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APPENDIX A: CALCULATION OF THE 
EXCHANGE INTEGRAL 

The three dimensional integral in Eq. (|10|l can be writ- 
ten as a combination of integrals of the form: 



/>27T /> 1 f 1 

Rn{o>) = da xdx / ydy 
Jo Jo Jo 



1 sign(n) — (— 1)™ cosa 



y/ x 2 + y 2 a 2 — 2xyacosa 
(Al) 

where: n = 0,1,2, sign(n) gives the sign of n and 
sign(0) = 0. The values of the functions R n (a), for a = 0, 
are R (0) = and i? x (0) = R 2 (0) = w. We also have: 



M 1 ) = ^(-2 + 7r(ln2 + l/2)+4C-7r(l + ln4)/2) 



~ 1.109, 
i?i(l) = 8/3 + i? (l) ~ 3.776, 

where C ~ 0.915966 is the Catalan constant. 
Assuming that < a < 1 we define: 



(A2) 



R n (a) 



da[sign(n)l - (-1)™ cos a]K(a, a) , (A3) 



where K{a,a) is given by: 



K(a,a) = — 7t — (1 + a 3 ) + (1 + a 2 )\/ 1 + a 2 — 2a cos a — (1 + a 3 ) cosaln(l — cosa) — a 3 cos a In a 
3ar L 

+ cos ol In (a — cos a + \/\ + a 2 — 2a cos a) — a 3 cos a ln( 1 — a cos a + \/\ + a 2 — 2a cos a) (A4) 

I 



This expression allows us to obtain the expansions: 



Ro(a) ~ -[-alna+S Q (a)], (A5) 
R n (a) ~ ^[3+ {-l) n alna + S n (a)}, (A6) 



for n — 1 , 2 and 

So (a) = 



45 



2 In 2- - I a a 3 - 

6/ 80 1792 



175 
18432 



a , 



S n (a) = 



(A7) 



-I -lr'i21n2-i)a-^ 2 + (-l)"A a 3 



3 a 4 + (-iy 



45 



15 



64 

y ' 18432 



1792 1024 



7 0> 



(A8) 
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Note that R\(a) — #2 (a) = 2R (a) is always satisfied. 



APPENDIX B: CALCULATION OF THE 
EXCHANGE ENERGY IN THE PRESENCE OF 
DISORDER. 

We write the one-electron energies in the absence of 
interactions and disorder as: 



The occupancy of a given state at fixed chemical po- 
tential, ep, is: 

n± = f F -lmG ± (k,oj)dLu, (B3) 

where a frequency cutoff, uj c is also defined. 

The total number o electrons, n, and the kinetic energy 
can be written as: 



e k = ±v F (k) , 



(Bl) 



up to some cutoff k c , where the two signs correspond 
to the two bands in the electronic spectrum. Using the 
CPA, the one electron Green's function can be written 
as: 



G ± (k,uj) = 



1 



w - E(w) - e± ■ 



(B2) 



ae\k\n" k \kdk . 



(B4) 



These one-dimensional integrals are calculated numeri- 
cally. Finally, the exchange energy is: 



£cxch = -^4 / d 2 fci / d z k 2 



2j [n+(fci) + n-(k 2 )} 2 + [n+(fei) - n-(fc 2 )] 2 cos[0(fc 1 ) - 0(fc 2 ) 



Ifci -fe 2 | 



(B5) 



r 



and: 



cj)(k) — arctan ( 



(B6) 



This expression can be reduced to a three-dimensional 
integral, which is calculated numerically. 

The total energy, E tot (n) = K(n) + E exch (n), can be 
written as: 



E tot (n) = £ to t(n T ) + Etot^) . 



(B7) 



The exchange instability towards ferromagnetism implies 
that: 

E tot (n/2 - 5n) + £ to t (n/2 + <5n) < 2S tot (n/2) , (B8) 
so that: 



9n 2 



< 0. 



n/2 



(B9) 
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